DFT-based many-body analysis of electron transport through molecules 
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We present a method which uses density functional theory (DFT) to treat transport through a 
single molecule connected to two conducting leads for the weak and intermediate coupling. This case 
is not accessible to standard non-equilibrium Green's function (NEGF) calculations. Our method is 
based on a mapping of the Hamiltonian on the molecule to a limited set of many-body eigenstates. 
This generates a many-body Hamiltonian with parameters obtained from ground state L(S)DA- 
DFT calculations. We then calculate the transport using many-body Green's function theory. We 
compare our results with existing density matrix renormalization group (DMRG) calculations for 
spinless and for spin-1/2 fermion chains and find good agreement. 
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I. INTRODUCTION 

The rapid development in the field of electrical trans- 
port through molecules and quantum dots has induced 
a considerable effort to investigate the physical mecha- 
nisms behind it. For a conceptual understanding of these 
phenomena it is indispensable to develop methods in- 
volving a minimal number of approximations. Different 
schemes have been used for the calculation of the con- 
ductance of such systems. The most popular ones are ab 
initio methods based on density functional theory (DFT) 
[1] in connection with the non-equilibrium Green's func- 
tion (NEGF) formalism which has been successfully used 
for understanding coherent transport through molecules 
in the strong coupling or off- resonant regime [2-4]. Also 
DFT is known to yield good results for ground state cal- 
culations. Although it has been argued by Stefanucci and 
Almbladh [5] that appropriate time dependent functional 
should also give good results for transport, the currently 
available functionals are not satisfactory for transport 
in the weak coupling regime where the Coulomb inter- 
action between the electrons dominates their dynamics 
[6]. Using the available ground states DFT functionals 
in particular gives a poor description of ionization, ad- 
dition and excitation energies and these states play an 
important role in transport. 

To illustrate the failure of DFT, we consider the usual 
case of non-ferromagnetic leads which always yields a so- 
lution in which spin up and down have the same occu- 
pation. Transport through a single level through which 
one or two electrons can flow, is described by a single- 
level Anderson type model. The most general form of the 
density matrix (restricted to the transport level) during 
transport, is given by: 

p = |0)(0| + a\ t)(t | +a|;><4 1 +6| U)(U I (1) 

Local density approximation (LDA) in DFT yields for 
this case a restricted solution, which does not distinguish 
between the last three terms and it is well known it can- 
not produce the correct step-like behaviour of the cur- 
rent as a function of voltage [7]. Instead, for each level, 
the current rises gradually with bias up to a maximum 



value. Part of the shortcoming of DFT(LDA) may be 
corrected for by adding a self-interaction correction (SIC) 
which first leads to a plateau corresponding to a single 
conduction channel before it steps to a next plateau at 
maximum current corresponding to the two conducting 
channels [8, 9]. Although the SIC method is an improve- 
ment over the standard DFT, it still gives wrong results 
for the current value through a single plateau: the SIC 
method predicts the current value of the plateau to be 
half of the maximum current while, as we will show in 
Sec III, it should be 2/3 of the maximum current. 
These shortcomings have induced the development of dif- 
ferent methods for weak coupling regime. As an ex- 
ample, the many-body effects that are not captured by 
DFT-NEGF can be obtained by the GW approximation 
method, which however is very time-consuming [10]. 
Combining DFT with rate equations can be used to de- 
scribe the electron transport in the weak coupling regime 
[11, 12] but this technique requires fit parameters and 
cannot show the broadening of the isolated levels due to 
the coupling (except for the temperature broadening). 
However, DFT is a powerful means to calculate the total 
ground state energies and this leads us to exploit this 
advantage of DFT in this regime. Thus our purpose 
is to present a technique relying on the combination of 
DFT and many-body NEGF approach which deals with 
transport in the weak coupling regime. Our method com- 
bines local spin density approximation (LSDA) for differ- 
ent numbers of electrons with many-body Green's func- 
tions (GF) to calculate the transport through a molecule, 
weakly connected to two non-interacting leads. We illus- 
trate our method using an interacting hopping chain for 
particles with and without spin. The latter case allows 
for a comparison with density matrix renormalization 
group (DMRG) calculations [13, 14]. We not only obtain 
excellent values for addition and ionization energies, but 
also good agreement of the location and the line shapes 
of these resonance levels when comparing with results 
based on DMRG method. The line shapes are the result 
of the coupling between the states on the molecule to the 
leads, which we also calculate using our DFT states. It 
is envisaged that the method of the paper will be useful 
within ab initio quantum chemistry calculations for elec- 
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tron transport. 

The organization of this paper is as follows. In section 
II the model for spinless and spin-1/2 fermions is defined 
and then our method is explained. The results for the 
single level inside or near the bias window are discussed 
in Sec III. Then the results for the more complicated 
case with two levels inside the bias window are presented 
in Sec IV. The conclusions in Sec V briefly summarize 
our ideas. The appendices include further details con- 
cerning Bethe-Ansatz solution for spinless fermions (A), 
L(S)DA-DFT for the Hubbard model (B) and calculating 
the transport through a Coulomb island (C). 
In this paper, we use the term 'level' to indicate a chem- 
ical potential corresponding to an energy resonance on 
the molecule. 



II. MODEL AND METHOD 
A. spinless fermions 

The systems studied here consist of a small region 
where Coulomb interactions are present, weakly coupled 
to two non-interacting, semi-infinite leads (see Fig. 1). 
The interacting region contains one or several quantum 
dots in series. The Hamiltonian of the entire system is 



H = H 



leads 



H, 



coupling 



H, 



molecule 



(2) 



The Hamiltonian for spinless fermions with interaction 
reads [15]: 

N L -i 

^molecule = —t ^ [d\di+\ + h.C.] 



i = l 



N L -1 



N L 



i=l i=l 

where rii = d\di and Nl is the length of the interacting 
chain. The parameter t represents the hopping rate and 
U describes the inter-site Coulomb interaction. The cre- 
ation and annihilation operators, d\ and di acting on site 
i satisfy the usual anticommutation relations. In addi- 
tion, the external gate potential, V g , can be applied to 
the interacting region which is included in the energy e. 
For the nonintcracting leads, 

N L -1 

-ffleads = - E tc E K^i+M + h.C.} (4) 
r)=L,R i=l 

where t c is the hopping term in the contact part, and 
the label r\ = L,R for left (L) and right (R) lead. The 



eigenstates are 



with energy [16] 



E = E — 2t c cos ka 



(5) 



U=0 



+V/2 




FIG. 1. A short Hubbard chain connected to two non- 
interacting leads. 



where 



e ika = _ Q 



q = -E- (6) 



We take a = 1. In addition, the bias voltage can be 

applied to the contacts. 

The coupling Hamiltonian reads 



H, 



coupling 



E tvL ; d j + h - c -\ 



(7) 



j'Gmolecule 



The Hamiltonian for the central part, i? mo iocuio, can be 
solved exactly using the Bcthc-Ansatz solution [17]. For 
such a system, Takahashi [18] gives the equations which 
should be solved for the density n. This is briefly ex- 
plained in Appendix A. 



B. spin-1/2 fermions 

Adding the spin as an extra degree of freedom, -£/i oa d s 
and -ffcoupiing are similar to the described Hamiltonians 
for the spinless case but a spin index a =f, ! is added to 
the creation and annihilation operators. 

Nl-I 

Pleads = - E tc E E l C i, V ,a C i+l,V,<? + h - C -\ ( 8 ) 
r]=L,R a i=l 

The Hamiltonian of the interacting region reads 

N L -1 

-ffmoloculc = — ^E E [^i,o-^+l,CT + h.c] 
N L N L 

+ u E d h d n4id ti + e E E 4a* (9) 



In this case, the Coulomb energy is on-site (as two parti- 
cles may now occupy the same site) and we parametrize 
it again by U. 
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FIG. 2. The states of the interacting chain in (a) are mapped 
onto a single dot which contains several interacting levels (b). 
The coupling between these levels and the leads are defined 
a.s t cH . The new model includes the intra-level Coulomb in- 
teraction and the inter-level interaction between the states. 



C. Method 

Our method for calculating fcrmion transport through 
the interacting chains starts by expressing the molecular 
Hamiltonian in terms of its (many-body) exact eigen- 
states \S): 

H molecu i e = J2\S)E s {S\ (10) 

s 

Because of the two-body character of the Coulomb po- 
tential, we can formulate this Hamiltonian in terms of 
creation and annihilation operators for (spin-) orbitals 
| a) with a Coulomb interaction: 

-^molecule = £ a (ft a d a + - ^ U a pdld a <tfpdp (11) 

Note that the quantum number a includes the spin (for 
spin- 1/2 particles). The eigenstates \S) are then the 
states \n a ), where n a = 0,1 represents the occupation 
of all (spin-) orbitals \a). Note that Eq (11) is a refor- 
mulation of the original Hamiltonian (Eq (3) or (9)) in 
terms of an interacting multi-level interacting Anderson 
model. Eq (10) is another reformulation of these Hamil- 
tonians. We shall use Eq (11) to find the specific form of 
the many-body eigenstates S appearing in Eq (10). 
Our method is based on a mapping of the Hamiltonian on 
the molecule to a limited set of many-body eigenstates 
(Fig. 2). We first find the set of parameters {e ai U a ^) 
from DFT ground state calculations. For the interacting 
chains used in this paper, we use the LDA parameteriza- 
tion based on the Bethe-Anstaz solution for interacting 
fermion chains [17-21]. In the case of spin-1/2 particles, 
we have used an accurate LSDA parameterization, given 
by Franga, Vieira and Capelle (FVC) [22]. 

We first consider particles without spin. DFT allows 
for calculating ground state energies for any number of 
particles. For a chain of N sites as in Fig. 2a, this num- 
ber varies between and TV, so DFT gives us TV + 1 



energies. Considering the system as in figure Fig. 2b, we 
would need N chemical potentials e Q and N(N — l)/2 
Coulomb interactions between these levels which adds 
up to N(N + 3)/2 parameters. It is therefore clear that 
this parametrization is highly non-unique. The situation 
is different in the case of spin-1/2 particles. For par- 
ticles with spin 1/2, we can vary the particle number 
M between and 2N and we can vary the polarization 
(A/ t - Mi) between -Min(M, N) and Min(M, N) yield- 
ing N(N + 3)/2 different ground state configurations pre- 
cisely the number of parameters needed in the Fig. 2b. 



For small bias voltage, the transport is dominated by 
a single chemical potential, corresponding to a transition 
from N to N + 1 particles. For this case, we can still 
apply our method to spinless particles. 
We calculate transport for the system consisting of a 
molecule described by the Hamiltonian (11), coupled to 
the noninteracting leads. It is therefore necessary to eval- 
uate the coupling for the different (spin) orbitals \a) to 
the leads. We do this by projecting the original chain 
Hamiltonian (3) or (9) onto two many-body states of 
the isolated central region, differing by one particle. We 
call those states |Sjv-i) and \Sn). Sn is obtained from 
Sjv-i by putting a particle into level a which is empty in 
Sjv— l, \Sn) = \Sn-i,&)- We explain the method for the 
spinless case. The Hamiltonian, formulated in the space 
spanned by \Sn) and |5jv-i) is 

H = \S N -i)E SN _ 1 (Sjv-i| + \Sn)Es N {Sn\- (12) 

The coupling Hamiltonian describes the hopping of a par- 
ticle from or onto the left lead to the leftmost site of the 
central region and a similar description can be used for 
the right lead. We anticipate that the effective coupling 
of a level a to the leads varies with the amplitude of 
that state on the left- and rightmost sites respectively. 
We calculate this effective coupling using the projection 
operator: 



P=\S N - 1 )(S N - 1 \ + \S N )(S N \ (13) 

The process in which a particle hops from the left lead 
onto the leftmost site of the central region is described 
by the following term of the coupling Hamiltonian (the 
calculation for the hopping to the right lead is similar): 

^coupling = t L d\c L . (14) 

The full Hamiltonian with the central region projected 
onto the subspace spanned by Sn-i and Sn, contains 
transitions of the form i? C oupling = t cff IS'jv) (Sn-i \cl- 
Projecting i/coupiing onto the span of Sn and Sn-i gives 

P^couplingP - t L (\S N )(S N \)d\c L (\S N ^)(SN-i\il5) 
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using c([\Sn-i) = |1; jSjv— i) leads to 

P -ffcoupling-P = i L|S'Ar)(S'Ar|l; SjV-l) vSjV-l \ c L = 

tf\S N )(S N -i\c L (16) 

which requires i^f = ti(5jv|l; Sjv-i)- Here L is the site 
of the left lead connected to the central region, and 
|1;Sjv-i) denotes an antisymmctrized state obtained by 
adding an electron on site 1 to a central region containing 
N — 1 particles in state |5jv-i)- 

In DFT, the approximated cigenstate is given as 

\Ss) = -±= y £T,p\<p$ l ...<f$ s ), (17) 
V-iV. p 

i.e. a Slater determinant composed of the single-particle 
DFT orbitals ip^ found within the N particle ground 

state (y ] is a sum over permutations and r\p is the sign 
p 

of the permutation). Defining l^ -1 ) = |1), t^ ff reduces 
to 

tf a = tL^lpU^t^pJ = *i x dct (S) ( 18 ) 

P n 

where S is the "overlap matrix", Ski = {^k ~~ Ivf ) and a 
denotes the highest orbital of Sjv- In the case of spin-1/2 
particles, the calculation of the effective coupling depends 
on tpt, ip l which leads to t cS = t L)R x det(5 t ) x det(5 ,J '). 
The effective coupling mainly depends on the shape of the 
orbitals, in particular their values on the outermost sites 
of the molecule. However, this shape for an electron with 
spin-up also depends on whether a spin-down electron 
occupies the level. We account for this by writing, for 
the coupling of a spin-up electron 

tf,R{n X ) = (1 - H)tfJ{ H = 0) + n it l%\ ni = 1) 

(19) 

It turns out that the values of the coupling for the 
two occupations = and = 1 differ only slightly 
(less than 2 percent). We neglect the influence of the 
occupation of the other orbitals. 

Our method for calculating the transport now consists of 
the following steps: (i) Calculate the ground states of the 
molecule for different charge states N and polarizations 
p = M t - Mi using L(S)DA-DFT. (ii) Infer the values 
for e a and U a p from these results. (iii) Calculate 
the effective coupling for the (spin-) orbitals \a). (iv) 
Calculate the transport for the Hamiltonian (11) coupled 
to non-interacting leads by an a— dependent coupling 
obtained in step (iii). 

The retarded and advanced GFs, G r and G a , for the 
transport calculation can be derived from the equation 
of motion. In order to find the lesser GF, G < , we use the 
Kadanoff-Baym equation 

Go 1 G< = £ r G< + £<G Q (20) 



For details see Appendix C. Once these GFs are known, 
the current can be calculated from a Landaucr type of 
equation 

I=jj Tr{ r ^-(G''-G Q )}(/( W , Aii )-/( W ,^))^ 

(21) 

where Tj — i(£J — £J ) and is the retarded self-energy 
and f(ui, Hj) is the Fermi distribution of lead j. 

A few remarks arc in order. In practice, we select only 
a limited set of many-body states, notably those whose 
charge additions and ionizations correspond to a chemi- 
cal potential inside or near the bias window. This means 
that we neglect the low-lying (spin-) orbitals (that are 
always occupied) and the higher orbitals (that are never 
occupied). This enables us to treat the spinless fermion 
transport with low bias, even though we cannot find all 
the values e a and U a p in that case. 
Calculating the transport is a standard problem for 
a single orbital inside the bias window for spinless 
fermions. However, approximations are necessary as soon 
as Coulomb interactions become relevant. We follow the 
simplest approach, in which correlation with the leads 
are neglected [23, 24]. This means that we will not ob- 
serve the Kondo resonance for spin-1/2 fermions. More 
elaborate schemes are possible, in particular slave-boson 
techniques which do take these correlations into account 
[25]. 

Even within our approach, transport through a single 
orbital for spin-1/2 fermions is already nontrivial. For 
more orbitals, several schemes based on further approxi- 
mations have been devised (see e.g. B. Song et al. [26]). 
We treat the full problem of spin-1/2 transport through 
two orbitals (four spin-orbitals) neglecting only the cor- 
relation with the leads. For details see Appendix C. This 
already allows for 8 transport channels (7 in the case of 
degenerate levels). In molecular electronics, bias voltages 
are hardly ever high enough to observe that many states, 
so we do not consider larger systems. 
A similar approach has been proposed by Yeganch et al. 
[27] based on quantum chemistry calculations for ground 
states and excited states. Also a time-dependent version 
of LDA functional for a similar model has been used by 
Kurth et al. [28] to investigate the transport within time- 
dependent DFT. Other approaches to describing trans- 
port in the weak coupling limit were based on the configu- 
ration interaction method for the central region, in com- 
bination with rate equations [29, 30] and with integra- 
tion over scattering states constructed through a Wigncr 
transform [31]. 



III. RESULTS FOR A SINGLE LEVEL INSIDE 
THE BIAS WINDOW 

In this section, we first present the results for the 
spinless fermions and compare with DMRG results 
obtained by other groups [15]. Then we discuss the 




FIG. 3. linear conductance of 7 noninteracting dots, U = 0. 
(a) our results and (b) the results from Ref [15] (Copyright 
Wiley- VCH Verlag GmbH and Co. KGaA. Reproduced with 
permission) .V b = 2 • 1CT 4 , t L>R = 0.5, t c = 1, t = 0.8. 
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FIG. 4. linear conductance of 7 interacting dots for weak 
(squares) and strong (circles) interaction, (a) our results and 
(b) the results from Ref [15] (Reproduced with permission). 

v b = 2 • nr 4 , t L ,R = 0.5, t c = i, t = o.8. 



result for spin- 1/2 particles. We consider small bias, 
so that at most one level lies inside or near the bias 
window. Energies and parameters with the dimension 
of the energy can from now on always be assumed to be 
given in Volts(V) and the current and conductance units 
are e/h and e 2 /h respectively. 



Mr 

28 +U 
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FIG. 5. two different situations to extract e a , U by FVC 
parameterization. 



B. Spin- 1/2 fermions 



A. Spinless fermions 

We first consider the results for U = 0. The linear con- 
ductance versus the gate voltage in the case of spinless 
fermions is shown in Fig. 3a for 7 non-interacting sites 
compared to the Fig. 3b obtained from Ref [15] based on 
the DMRG method. The gate voltage is applied in order 
to shift different resonant levels across the narrow bias 
window. For this case, the peak locations are easy to 
get at the right position, i.e. — 2<cos(fca), where t is the 
inter-dot coupling and k is the wave vector that fits on 
an isolated chain of 7 dots. The agreement between peak 
widths in the two figures shows that our way of calculat- 
ing the effective coupling seems correct. The last peak is 
higher in our simulation than in the DMRG result, due 
probably to a limited number of V g values used in the 
latter. The linear conductance versus the gate voltage 
for 7 interacting sites in the cases of weak (U/t = 1) and 
strong (U/t = 3) interactions are shown in Fig. 4a. Since 
this is also in agreement with DMRG results in Fig. 4b, 
we conclude that our method is reliable. In both figures 3 
and 4, applying the negative gate voltage, will show three 
peaks for the linear conductance in that region. 



For spin-1/2 particles, the ground state energies for 
one site containing one and two electrons have been cal- 
culated using the FVC parameterization which gives the 
values for e and U (Fig. 5). The many-body approach de- 
scribed here can then be applied to calculate the current 
through one or several quantum dots. The result for one 
dot is shown in Fig. 6 for Ep — 0, e — 0.5 and U = 0.2. 
Our hybrid method gives two steps at the expected posi- 
tions V = 2|e - E F \ = 1.0 and V = 2\e - E F + U\ = 1.4. 
We compare this calculation with the LDA-DFT-NEGF 
method, using the LDA parameterization by Capelle et 
al. [20, 21] (see Appendix B). The LDA curve gradually 
increases from V — 2|e — Ep\ till about 1.21^ where it 
reaches a level corresponding to transport through both 
channels. We see that including the spin explicitly into 
the transport calculation makes a substantial difference, 
even though on average the spin on the molecule is zero. 
Therefore, using a restricted exchange correlation func- 
tional is bound to give wrong results. 

As we explained in the introduction, the SIC predicts 
the value of half of the maximum before it steps up to 
its maximum value [8] , while the correct value of the cur- 
rent in the weak coupling limit at this region is 2/3 of the 
maximum current. This can be found using rate equation 
calculations [32] and it can be understood from the fact 
that there are two one-electron channels (corresponding 
to spin up and down), but only one two-electron chan- 
nel. 

The negative slope after the second step is due to the 
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FIG. 6. current through one quantum dot for LDA and for 
the many-body combined with LSDA (FVC parameteriza- 
tion) for nn,L = E F ± V/2, E F = 0, e = 0.5, t L ,fl = 0.1, 
t = 1.0 and U = 0.2. 




FIG. 7. Conductance of three quantum dots system compared 
to the results of embedded-cluster approximation (ECA) of 
Ref [33]. £7 = 1, t LM = 0.3, t c = 1, t = 1, V Mas = 0.006. 



fact that the density of states in the leads is not con- 
stant. Therefore, at different biases, the leads supply a 
different number of electrons. Indeed, in the case of wide 
band limit (where self-energies are independent of energy 
and bias voltage) the negative slope disappears. 
We have also compared our results for three coupled dots 
with DMRG in combination with the embedded-cluster 
approximation (ECA) of Ref [33] in Fig. 7. The resonance 
peak position is precisely in agreement with the Fig. 12c 
of Ref [33] and also the general line-shapes of the peaks 
are very similar to that where two central peaks are wider 
and four peaks at two sides are narrower than the central 
ones but the lowest values between peaks are different. 
Fig. 12c of Ref [33] is made assuming even numbers of 
electrons in the contacts and hence neglects the correla- 
tions held responsible for the Kondo resonance. Although 
difference with Fig. 12a of Ref [33], which is believed to 
be the correct, seems significant, we expect them to be 
much less pronounced at higher bias. As it is, the agree- 
ment shows that within the approximation made in our 
Green's function approach, our method gives the correct 
prediction. An improvement which would include these 
correlation will be considered in future work. 
Fig 8 shows the occupation of the Hubbard site (n^+ni) 
versus applied bias voltage for five coupled dots with 
Ep = 2.4. We find e = 2.54 and U = 0.18 and we see the 
steps at V = 2|e - E F \ = 0.14 and V = 2\e - E F + U\ = 
1.1. The density of the level at zero bias is non-zero which 
is due to the broadening of lowest unoccupied molecular 
orbital (LUMO) which contributes in the transport. We 
have shown the differential conductance M? as well for 

aV 

five dots which can be compared with DMRG results for 
spinless case [15]. 




FIG. 8. Left: Five coupled dots system. E F - 2.4, U = 1.0, 
ih,R = 0.4, t c = 2, t = 1.0. Extracted values are e a = 2.54, 
£7 = 0.18, tf R = 0.1, tf R = 0.1. Right: 

differential conductance for five quantum dots. E F — e a — 
2.5479, U = 1.0, t LR = 0.5, t c = 1. Extracted values for 
coupling are tf^ N+1 = tf tR _ N+1 ^ N+2 = 0.143375 and £7 = 
0.184. The shape of the curve is in agreement with DMRG 
results (Ref [15]). 



IV. RESULTS FOR TWO LEVELS INSIDE THE 
BIAS WINDOW 

For most experimentally relevant situations, the 
single level problem with interaction will be adequate. 
However, if the molecule possesses a symmetry (which is 
not destroyed by an imbalance in the contact geometry) 
molecular orbitals may become degenerate. Indeed, for 
chains with more than one site, we find degeneracies in 
the spectrum of the isolated molecule. 
In this section we therefore consider a problem with two 
degenerate levels (see Fig. 9) which may lie inside the 
bias window. For this case, we do not know of reliable 
calculations to compare our results with. However, in 
view of the good agreement with DMRG for the single 
(interacting) orbitals, we expect the results presented 
here to be reliable. We label the four states for two 
levels 1-j-, lj, , 2-f, 2^ respectively. We map our system 
(shown in Fig. 1) to a model with two orbitals with 
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FIG. 9. Schematical model for a two level system. 
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FIG. 10. Five proposed configurations to extract the values 
of ei, £2, Un, U22 and C/12. 



chemical potential that may be degenerate and include 
the intra-level Coulomb interaction (Un and U22) and 
the inter-level Coulomb interaction (Uw). We assume 
that the inter-level interaction does not depend on spin 
(Fig. 9). In this case we should calculate £1, £2, Uxi, 
U22 and U\2- For this purpose, we calculate the ground 
state energy by FVC parameterization in the five cases 
shown in Fig. 10. Thus by having these energy values 
we can calculate the mentioned energy levels and the 
Coulomb interactions. From these values, we can then 
investigate the transport (sec Appendix C for details 
about the method). 

An important feature of our method is that it can 
produce I-V characteristic for rather high bias voltages. 

We start with a chain consisting of two dots. As 
explained in Sec. II C, we first map the spectrum of this 
chain onto two one-electron levels that can be occupied 
by spin- up and/or spin-down electrons (see Fig. 9). 
This is done by calculating the ground state energies 
for all the configurations shown in Fig. 10. From these, 
we find the appropriate e and U values. The results 
are shown in Fig. 11 for a chain with parameters given 
in the caption of that figure. The transitions seen in 
the curve of Fig. 11 are displayed in Fig. 12, together 
with the predicted energies for these steps based on 
the parameters of the two- level model of Fig. 9. The 
first and second steps (a) and (b), will take place when 
the bias voltage is not high enough to encompass both 
levels but it is high enough to cross the first level. As 
this level can be occupied by two electrons we see two 
steps between V& = and Vb = 3. A higher bias voltage 
enables the occupation of the second level. The third 
step shows the addition of an electron to the second 
level. 

We also have mapped a Hubbard chain of three 
interacting dots onto the 3-lcvcl model. To extract 



FIG. 11. Occupation of a two quantum dot chain with t = 1, 
t c — 6, th,R = 0.3, V g — 1.4, intial U = 1 which lead to 
ei = 0.4, e 2 = 1.86, Uu, 22 = 0.46, U 12 = 0.53. tf 



it->H 



= t 



H- 



_ j-eff 
*2t — I 2t- 



eff 
l 0->lt 



2 j_ = 0.212132. The red curve shows 
the occupation of the first level while the green one shows the 
occupation of the second level. The blue one is the sum of the 
occupations. The six steps shown in this curve correspond to 
different transfer process presented in Fig. 12. 











\ r 




4.18->7 > 5 6 

ti / V=2(3.38|=6.76 \ ^ 



FIG. 12. Different transfer process of electrons corresponding 
to six consecutive steps shown in Fig. 11. (a) shows the trans- 
port of the first step in density curve (b) shows the second 
step and so on. Other processes play a role but they do not 
show up as separate steps in this graph. 



the nine values of e 1: e 2 , £3, ^11,22,33 and ^12,23,13; we 
considered the nine ground state configurations shown in 
Fig. 13. Here we take the two lowest levels to be inside 
or near the bias window. The result is shown in Fig. 14. 

We have implemented our method for at most two 
levels inside or near the bias window. However, it is 
possible to use the method for more than two levels 
inside the bias window which makes the computation 
time-consuming due to the larger dimension of matrices. 
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V. CONCLUSIONS 

In conclusion, we have proposed a method based on 
DFT which can accurately predict the transport in the 
weak coupling regime through an interacting chain. In 
our approach, we map the interacting part of the system 
to several interacting energy levels and take the Coulomb 
interactions into account. The dot occupations show dif- 
ferent steps corresponding to different transfer processes 
of electrons from the leads to the interacting region. Our 
method is a new opening for using DFT first principle cal- 
culations to investigate the transport through molecules 
in the weak and intermediate coupling limit. We do not 
have the observation of Kondo in our method but it can 
be included using more advanced approximations. We 
plan to implement our method into a quantum chemi- 
cal DFT code to calculate transport through experimen- 
tally relevant devices. As ground state DFT can predict 
excitation energies reasonably well (provided a separate 
self-consistent calculation is performed for each particle 
number and polarization) we should be able to reveal 
several energy levels. Furthermore the good results for 
the peak broadcnings in our model system arc promising, 
although the coupling strengths in experimental devices 
suffer from sample to sample variation. The same holds 
for the dielectric environment which may affect the lo- 
cation of the levels substantially. For this, and for the 
alignment of the levels to the Fermi energies of the con- 
tacts, a constrained DFT approach may be useful. 
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FIG. 14. Occupation of a three quantum dot chain with Vg = 
1.6, U = 1.2, t L , B . = 0.4, t c = 6, t = 1. Extracted effective 
couplings are tf^^^^ = 0.2, tf, R = 0.275279, 



tf R = 0.282842. 



equations are 



and 



where 



p(k) dk 

-B 



l = 27rp(*)- / T(k,q)p(q)dq 

-B 



T(k,q) = —6(k,q) 



and 

8(k,q) = 2tan~ 1 
E is then given as 



gsin(V) 



cos(^)-|cos(^) 



E=---2tjjoo sk+ -)dk 



(Al) 



(A2) 



(A3) 



(A4) 



(A5) 



Appendix A: Bethe-Ansatz solution for spinless 
fermions 

Here wc briefly describe the numerical approach for 
finding the exact ground state energy. Takahashi [18] 
gives the equations which should be solved for the density 
n, the 'quasi-momenta' k and the integer limit B. These 



and the exchangc-corrolation energy is then defined as 

2 



E xc = E- E{U = 0)-U[n 



1 



(A6) 



Since the function 8 depends on the momenta, the prob- 
lem has to be solved self-consistcntly. 
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FIG. 13. Nine proposed situations to extract the values of e\, 

£2, £3, £Al,22,33 , U12, 23,13- 



Appendix B: L(S)DA-DFT for the Hubbard model 

In order to construct a DFT Hamiltonian with param- 
eteres U, t for a Hubbard chain based on L(S)DA, an ex- 
pression for the exchange-correlation potential is needed. 
This potential is based on the exact ground state energy. 
The exact ground state energy e(n,m,t,U) (for density 
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n = + ni and magnetization m = — nj,) of the 
Hubbard model can be obtained using the Bethe-Ansatz 
[19, 34, 35]. An approximate analytical expression for the 
unpolarized case (m = 0) was proposed by K. Capelle et 
al. [20, 21]. This reads 



e(n < l,m = 0,t,U) 



2t(3(U/t) 



sin 



p{U/t) 



(Bl) 



where n = N/L is the Hubbard site occupation, TV, L are 
the number of electrons and Hubbard sites respectively 
and (3 is a function of the ratio U/t. It can be determined 
from the implicit equation 



2tgU/t) . 
sm( 



p{u/ty 



-At 



Jq(x)Ji(x) 



(B2) 

here Ji=o,i{x) are the Bessel functions of the first kind 
[20, 21]. For n > 1, the energy is found from the particle- 
hole symmetry 

e(n > 1, m = 0, t, U) = e(2 - n, m = 0, t, J7) + U(n - 1) 

(B3) 

From the energy, we can obtain an analytical expression 
for the exchange-correlation potential which, in the un- 
polarized case, is 

V„(n,m = 0,t,U) = !%* = 
£ [e(n, m = 0, t, U) - e(n, m = 0, t, 0) - e ff (n, C/)QB4) 



where the Hartree-energy is 

e H (n,U) = Un 2 /4 



(B5) 



In the case of non-zero magnetization, the energy expres- 
sion e(n,m,i, U) has been constructed by V. Franga, D. 
Vieira and K. Capelle [22]. 

The linear Hamiltonian matrix dimension for the polar- 
ized case with polarization M, is ((jv+m)/2) ((jv-M)/2)' 
while in the LDA case has the dimension L: 



( 2-Tll -r 



H 















\ 
J 



(B6) 

This LDA-based Hamiltonian replaces the actual poten- 
tial felt by an electron when it enters into a site by a time 
average of electron occupation at that site. In the LSDA 
Hamiltonian, this potential is a function of the two spin 
densities (nf, njj. 



Appendix C: Calculating the transport through a 
Coulomb island 



1. Single level inside the bias window 

We start with one level inside the bias window and 
we explain how the GF for the spinless case and for the 
case with spin can be derived from the equation of mo- 
tion (EOM). The time derivative of the d operator (for 
molecule) and the c operator (for contacts) are (see [23]) 

2C a — £ a CL a 

+ X] U al3d a np + ^ KikaCvkffc ( C1 ) 



r)=L/R 

q 



(C2) 



Here, a denotes the spin-orbital a and a a is the spin for 
this a. k labels the traveling wave states in the leads 
r/ = L, R where L and R stand for left and right. If wc 
consider only one orbital and neglect the spin, the term 
with U a /3 drops out of the problem. In order to find the 
current, we use non-equilibrium GF theory, which focuses 
on the one-particle GF on the molecule, defined as 

G af3 = -i(T{d a (t)dl(t')}) (C3) 

where T is the time-ordering operator 

T{A(t)B(t')} = 9(t - t')A(t)B(t') t 0(t' - t)B{t')A{t) 

(C4) 

T always moves the operators with earlier time argument 
to the right. 

After Fourier transformation of the time domain, taking 
the time ordering carefully into account, the following 
equation for the GF is found [23]: 

(Lj-e a )G a0 (Lj) = 6 a p+J2 t n^k^) ( C5 ) 

T}k 



where 



rS(t-0 = -* , <rc, toa (t)4(0> 



t U'\ 



(C6) 



Using the EOM for c n ka a (t), an equation for is found: 
(w - e vk )T^(t - t') = t„G aP {u) (C7) 
Using the second equation to eliminate T^, we arrive at 

(u - e a - Y, (uj))G a p = Sap (C8) 

where 



(C9) 
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is the sclf-cncrgy. The sclf-cncrgy has a real (Hcrmitian) 
part which has the effect of shifting the resonant energies 
and which reflects asymmetries of the densities of states 
near those resonances. The imaginary (non-Hermitian) 
part broadens the resonances, reflecting the hybridization 
of the states of the central region with those of the leads. 
Taking a = j3 and writing Sg(o;) = A(w) + iT/2, we can 
write 



1 



e a — u) — A — iY /2 



(CIO) 



which is a Lorentzian function. 

Next, we consider an interacting dot for spin-1/2 
fcrmions. In that case, the EOM leads to the following 
equation for the GF: 

{u-e a )G aP (u) = 6 a0 + J2 U ai G%{w) + Y,t n T a £{u) 

7^a r/k 

(en) 

while the equation for remains the same. We have 

(2) (2) 

introduced a new GF G a ^, which is defined as G a Jp = 

-i(T{d a (t) ni {t)dl{t r )}) with n 7 (t) = dj(t)d\(t). This 
GF satisfies an EOM: 

(u-e a - U a -f)G^p = {n-f)5 a p 

+ £(*; r Sf + ^f-t;r^f) (ci2) 

r/k 



where 



rSf = -i{T{c nkaa (t)n 7 (t)d}(t')}) 



(C13) 



ifjf = ^(T{c nkay (t)d a (t)d 7 (t)dl(t>)}) (C14) 



t w\ 



rgf = -*(T{c nk ^(t)d\(t)d a (t)dl(t>)}) (C15) 

We now neglect correlation between the central region 
and the leads by keeping only T^^f which we approx- 



imate as r 



(2)a<3 
l,r}k 



Note that this mean field 



approximation only concerns the coupling between the 
central region and the leads, but not the Coulomb cor- 
relations within the central region. Thus by substituting 
G<® from (C12) to (Cll) and eliminating T by an equa- 
tion like (C7), it yields 



G aa (u}) — 

u -e a - (1 - (np))U 



(ui-e a - U)(uj - e a ) - S r [w - e a - (1 - {np)) 
where 57 = + T, r R and 



16) 



E 3 » 



-z (uj) 



(C17) 



and Im(^) > 0, Zj = -qj ± ^q) - 1, qj = E £ V/2 

and Ep is the Fermi energy of the grounded lead. 

To calculate the density sclf-consistently, the calculation 

of the lesser GF is also required 



/ Tp: dw 



(C18) 



which can be found from the Kcldysh equation Eq (C19). 

G< Q M = G r aa {u)^<{Lo)G a aa {u) (C19) 
G a is the advanced GF and the lesser self-energy is 



j=L,R 



(4i c 



(C20) 



Once the retarded and advanced GF are known, the cur- 
rent can be calculated from a Landauer type of equation 

I=jj Tr{ r ^-(G''-G a )}(/( W , AiL )-/( W ,^))^ 

(C21) 

where Tj = i(SJ - Ttf). 



2. Two levels inside the bias window 

Here we explain the transport calculation for the case 
we have two levels inside the bias window. Once we 
have calculated the energy values and Coulomb interac- 
tion by FVC parameterization, we can use a many-body 
approach which is a generalization of the technique de- 
scribed above. For the one-particle GF, we find the same 
equation as above: 

(u - E a )G a p = 5 al3 + ^2 U a~f G a7/3 + 



^ ' ^aa'G a 'p. 
a' 

where £ aQ < is the self-energy 



ijk 



(C22) 



(C23) 



For G (2) we obtain the EOM 

(ui — e a - U ay )G^p = (n^)5 a /3 
+ E U ^G%p + E Z aQ ,G^ P (C24) 

(5^a,7 a' 

Eq (C24) is not closed as a new GF, G a ^ S p, is generated 

(2) (3) 

in deriving the equation for G a ^p. The new GF, G a ^ S p, 
is 



G% p = -i(T{d a (t)n 7 (t)n s (t)d}(t')}) (C25) 



t /V\ 
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The EOM for G^gp reads 



(u — e a - U ai - U aS )G K a ^ S p = (n 7 na)(5 aj 9 

+ E U »«G$W + E ( C26 ) 



€7^0:, 7 a' 

which introduces another new GF, 



G«Lb = -i(T{d a (t)n 1 (t)n s (t)n e (t)dUt')}) (C27) 



for which the EOM is 



(w — £ - C/ Q7 - J7 a< 5 - J7 Q , e )G^ (5c(3 = (n 7 n$n e )5 a p 



EE aa >Gg 4 , 



(-1) 



(C28) 



We now must solve the set of equations (C22), (C24), 

(C26) and (C28) for the GFs G Q(3 to G^ Se0 . Wc organise 
these GFs into a 340 x 4 array 



G\p — {G a p,G a ,' ll3 ,G y a u 1 , S p 1 G y a , > „ 1 „ slef} ) 



(3) 



,(4) 



(C29) 



As all indices a, /3, ... run over four states, it is easy to see 
that the first index of this array runs over 4 + 16 + 64 + 
256 = 340 values. The equation for Q can be written in 
the form 



G^Q = (n) + Eg 



(C30) 



Here, C/ 1 is a 340 x 340 matrix, which, in the frequency 
domain assume the form 



At this stage, the process of generating new GF stops, as 
the EOM for G^ does not generate higher-order GFs. 



UJ-E a > - U a > 



oj — e a " ! — Ua'/'Y' — U a '"S' — U a f> 



(C31) 



and (h) is a 340 x 4 array 

(5 a p, (n-y}8 a '/3, (n-f'ns}5 a »i3, (n 7 »n, 5 ,n e )(5 Q "^) T (C32) 

and E is the 340 x 340 array with elements Eq, a . Qa , , 
where a\ denotes the index a of the composed index 
A = (a,a'j,a"j'5 7 a'"j"5'e). 

In order to find the lesser GF Q < , from which (n) can 
be found, we should use a Kcldysh or Kadanoff-Baym 
equation. These equations are conveniently derived from 
the Langreth rules [36]. These rules apply to the GF Q 
which is found from Eq (C30). Therefore, using the no- 
tation of that equation, the Kadanoff-Baym equation can 
be written as 

g^ l g< = ^ r g< + ^<g a . (C33) 

where g a is found as 

g a = (g^ 1 -et 1 ^). (C34) 



In electron transport theory, the Keldysh equation, 

g< = ((n) +g r z r )g<((n) + g a z a ) + g r ^<g a . (C35) 

is often used, with only the last term on the right hand 
side, as it can be shown for transport through a single 
channel, the first term vanishes for the single particle 
GF G Q( g. However, this is not the case when the 'higher' 
GFs G^ etc. are included (this was also pointed out 
by Song et al. [26]). For the calculation of the integra- 
tion in Eq (C18) one has to calculate the inverse of Go, 
many times (depending on the number of the integration 
points and the number of the required iterations to solve 
the problem self-consistcntly) , making the computation 
time-consuming. Therefore one could think of using the 
following approximations which cause reduction of the 
matrix dimension to 84 x 84 and 20 x 20 respectively. 

G% t[3 = -i{T{d a (t)n 7 (t)n s (t)n e (t)4(t')}) s 

ll(^)G { l p + (ns)G% + (n e )Gi% p ] (C36) 
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FIG. 15. Occupation of a two quantum dot chain. The 
chosen parameters are ei = 0.4, E2 = 1.8, Umtra-ievei = 
Uinter-ievei = 0.5, Y l ,r = 0.05, using wide band limit. 



G% S p = -i(T{d a (t) ni (t)n 5 (t)dl(t')}) ~ 

\[(n 5 )G% + { ni )G { X] (C37) 

Fig. 15 shows the effect of the approximations on the 
occupation. The curve corresponding to dimension 20 x 
20 shows four steps in agreement with the full GF while 
the one based on dimension 84 x 84 depicts five correct 
steps and finally six steps has been gained from the exact 
solution (dimension 340 x 340) as we discussed. However, 
for the lower biases the results based on approximations 
are still valid. 

As we can see in Fig. 15 the density does not exceed 2 
while in Fig. 11 the occupation exceeds 2 and this can be 
explained by the difference between the self-energies used 
in these two figures. The wide band limit has been used 
in Fig. 15 which supplies constant self-energies while in 
Fig. 11 the self-energy reflects the density of states in the 
leads not being constant. Therefore, at different biases, 
the leads supply a different number of electrons. 
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